library(tidyverse)
library(survival)
library(survminer)
library(readxl)

plot_my_surv2 <- function(fit, title){
  ggsurvplot(fit = fit, pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = title)
}

# blood 2020 Supp -------
blood.sup <- read_delim('mission/DLBCL-Btk/data/blood-supp.txt') |>
  select(contains(c('pid','btk','myd88','cd79b','marcks')))

blood.clin <- read_delim('DLBCL-Btk/data/blood-clinic.txt') |>
  left_join(blood.sup)

blood.var <- read_xlsx('DLBCL-Btk/data/blood2020-suppl2.xlsx', sheet = 2) |>
  filter(Gene %in% c('BTK','MYD88','CD79B')) 

blood.var |>
  count(Gene, p.annot) |>
  filter(Gene == 'MYD88') |>
  mutate(p.annot = case_when(p.annot == 'MYD88:p.*205R' ~ 'MYD88:p.L265P',
                             .default = p.annot),
         p.annot = fct_reorder(p.annot, n, .desc = TRUE)) |>
  ggplot(aes('', n, fill = p.annot, label = p.annot)) + geom_col() +
  coord_polar(theta = 'y', clip = 'off') +
  scale_y_continuous(labels = NULL) +
  scale_fill_viridis_d(option = 'turbo', begin = .1, direction = -1) +
  theme_pubr(legend = 'right') +
  labs(fill = 'CD79B mutation in DLBCL', x = '', y = '')

bld2020 <- blood.var |>
  filter(str_detect(p.annot, 'MYD88:p.*205R|CD79B:p.Y197[A-Z]')) |>
  mutate(p.annot = case_when(p.annot == 'MYD88:p.*205R' ~ 'MYD88:p.L265P',
                             .default = p.annot) |>
           str_replace(':p\\.', '_'),
         value = 1) |>
  pivot_wider(id_cols = `Patient ID`, names_from = p.annot, values_from = value, values_fill = 0) |>
  right_join(blood.clin, join_by(`Patient ID` == PID)) |>
  mutate(across(2:6, \(x)ifelse(is.na(x), 0, x)))

blood.clin |>
  survfit(Surv(OS_time, OS_status) ~ BTK, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with BTK mutation (Blood, 2020)")

blood.clin |>
  survfit(Surv(OS_time, OS_status) ~ CD79B, data = _) |>
  ggsurvplot(pval = TRUE,
             palette = c('blue','red'),
             legend.labs = c('WT', 'Mutation'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with CD79B-Y197H mutation (Blood, 2020)")

blood.clin |>
  survfit(Surv(OS_time, OS_status) ~ MYD88_265, data = _) |>
  ggsurvplot(pval = TRUE,
             palette = c('blue','red'),
             legend.labs = c('WT', 'Mutation'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with MYD88-L265P mutation (Blood, 2020)")

blood.clin |>
  survfit(Surv(OS_time, OS_status) ~ CD79B+MYD88_265, data = _) |>
  ggsurvplot(pval = TRUE,
             palette = c('blue','orange','green2','red'),
             legend.labs = c('WT', 'Mut-MYD88', 'Mut-CD79B', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients OS with MYD88-L265P & CD79B-Y197H mutation (Blood, 2020)",
                              width = 60))

bld2020 |>
  select(-c(age_gt_60:cluster_ICL)) |>
  write_csv('DLBCL-Btk/results/blood2020-survival.csv')

# Cell 2017 -------
gene_pan <- read_xlsx('mission/DLBCL-Btk/data/cell2017-supp.xlsx', sheet = 6, skip = 2) |>
  select(-c(2:3)) |>
  filter(Genes %in% c('MARCKS')) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'mutation') |>
  pivot_wider(names_from = Genes, values_from = mutation)

clin_meta <- read_xlsx('DLBCL-Btk/data/cell2017-supp.xlsx', sheet = 1, skip = 3) |>
  select(`Sample  ID`,`Overall Survival years`,Censored,`ABC GCB (RNAseq)`,`age at diagnosis`,Gender)

gene_var <- read_xlsx('DLBCL-Btk/data/cell2017-supp.xlsx', sheet = 4, skip = 3) |>
  filter(Gene.refGene %in% c('MYD88', 'BTK', 'CD79B')) |>
  select(Gene.refGene, AAChange.refGene, `648`:last_col())

our_var <- gene_var |>
  filter(str_detect(AAChange.refGene, 'L265P|Y197|Y196')) |>
  mutate(aa = str_extract(AAChange.refGene, 'L265P|Y197H|Y197S'),
         mutation = str_glue('{Gene.refGene}_{aa}', )) |>
  relocate(mutation)

pan_final <- clin_meta |>
  filter(!is.na(Censored) & !is.na(`Overall Survival years`)) |>
  mutate(sample = as.character(`Sample  ID`)) |>
  left_join(gene_pan)

var_final <- our_var |>
  select(mutation, where(is.numeric)) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'status') |>
  pivot_wider(names_from = mutation, values_from = status) |>
  right_join(pan_final)

pan_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ MYD88, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with MYD88 mutation (Cell, 2017)")

pan_final |>
  filter(`ABC GCB (RNAseq)` == 'GCB') |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ MYD88, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "GCB-DLBCL patients OS with MYD88 mutation (Cell, 2017)")

pan_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ CD79B, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with CD79B mutation (Cell, 2017)")

pan_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ BTK, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with BTK mutation (Cell, 2017)")

var_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ CD79B_Y197H, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with CD79B-Y197H (Cell, 2017)")

var_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ CD79B_Y197S, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with CD79B-Y197S (Cell, 2017)")

var_final |>
  survfit(Surv(`Overall Survival years`, !Censored) ~ MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with MYD88-L265P (Cell, 2017)")

write_csv(var_final, 'DLBCL-Btk/results/Cell2017-survival.csv')

# Nat Med 2018 -------
clin_meta <- read_xlsx('DLBCL-Btk/data/natmed-clinic.xlsx', skip = 1) |>
  select(individual_id, Gender, `Age-at first diagnosis`, PFS, PFS_STAT, OS, OS_STAT) |>
  filter(PFS != 'na' | OS != 'na') |>
  mutate(sample = str_replace_all(individual_id, '-', '_'),
         across(PFS:OS_STAT, as.numeric))

var_info <- read_xlsx('mission/DLBCL-Btk/data/NatMed-mutation-matrix.xlsx', skip = 1)

var_tidy <- var_info |>
  filter(Name %in% c('CD79B','MYD88','MARCKS')) |>
  select(-Description) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'status') |>
  mutate(status = ifelse(status >=2, 1, 0)) |>
  pivot_wider(names_from = Name, values_from = status)

var_final <- var_tidy |>
  mutate(sample = str_remove(sample, '_nullpair')) |>
  right_join(clin_meta)

var_final |>
  survfit(Surv(OS, OS_STAT) ~ MYD88, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with MYD88 mutation (Nat Med, 2018)")

var_final |>
  survfit(Surv(OS, OS_STAT) ~ CD79B, data = _) |>
  ggsurvplot(pval = TRUE,
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with CD79B mutation (Nat Med, 2018)")

write_csv(var_final, 'DLBCL-Btk/results/nat_med2018_survival.csv')

# mutation sum up ------
nat2018 <- read_csv('DLBCL-Btk/results/nat_med2018_survival.csv')

cell2017 <- read_csv('DLBCL-Btk/results/Cell2017-survival.csv')

bld2020 <- read_csv('DLBCL-Btk/results/blood2020-survival.csv')

nat_tidy <- nat2018 |>
  mutate(gender = Gender, age = `Age-at first diagnosis`, pfs_month = PFS, pfs_death = PFS_STAT,
         os_month = OS, os_death = OS_STAT, .keep = 'unused') |>
  select(-individual_id)

cell_tidy <- cell2017 |>
  mutate(os_month = `Overall Survival years`*12,
         os_death = as.numeric(!Censored),
         subtype = `ABC GCB (RNAseq)`,
         age = `age at diagnosis`,
         gender = ifelse(Gender == 'F', 'female', 'male'),
         sample = str_glue('cell2017_{sample}'),
         .keep = 'unused') |>
  select(-`Sample  ID`)

trini_tidy <- bld2020 |>
  mutate(sample = `Patient ID`, subtype = ifelse(cell_of_origin == 'UNC', 'Unclassified', cell_of_origin),
         os_month = OS_time/30, os_death = OS_status,
         pfs_month = PFS_time/30, pfs_death = PFS_status,
         MYD88_L265P = MYD88_265, MYD88 = as.numeric(MYD88_L265P|MYD88_noncan), .keep = 'unused') |>
  bind_rows(nat_tidy, cell_tidy)

trini_tidy |>
  relocate(sample:pfs_death, gender, age) |>
  mutate(CD79B_Y197X = as.numeric(CD79B_Y197F|CD79B_Y197S|CD79B_Y197H|CD79B_Y197N)) |>
  write_csv('DLBCL-Btk/results/sum3paper_survival.csv')

trini_tidy <- read_csv('DLBCL-Btk/results/sum3paper_survival.csv')

trini_tidy |>
  count(os_death)

## multivar cox -------
plot_forest <- function(data, type = 'OS') {
  data |>
    ggplot(aes(term, estimate, ymin = conf.low, ymax = conf.high)) +
    geom_pointrange() +
    geom_text(aes(label = str_glue('p={p.value}')), nudge_x = 0.2) +
    geom_hline(yintercept = 0, linetype = 'dashed') +
    coord_flip() +
    theme_pubr() +
    labs(title = str_glue('Multivariate Cox survival analysis on {type}'),
         x = 'Covariate',
         y = 'Hazard Rate (HR)')
}

coxph(Surv(os_month, os_death) ~ MYD88_L265P + gender, data = trini_tidy) |>
  broom::tidy(conf.int = TRUE) |>
  mutate(p.value = signif(p.value, 2),
         term = case_match(term, 'gendermale' ~ 'gender:male', .default = term),
         term = fct_relevel(term, 'age', 'gender:male')) |>
  plot_forest()

g1 <- last_plot()

## gross MYD88 ----
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ MYD88, data = _) |>
  plot_my_surv2("DLBCL patients OS with MYD88 mutation (n = 2153)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ MYD88, data = _) |>
  plot_my_surv2("DLBCL patients PFS with MYD88 mutation (n = 988)")

## MYD88-L265P with subtype facets ----
trini_tidy |>
  as.data.frame() |>
  survfit(Surv(os_month, os_death) ~ MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             data = as.data.frame(trini_tidy),
             facet.by = 'subtype',
             legend.labs = c('WT', 'Mutation'),
             palette = c('blue','red'),
             risk.table = 'nrisk_cumcensor',
             title = "DLBCL patients OS with MYD88-L265P mutation (n = 2153)")

## MYD88-L265P
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ MYD88_L265P, data = _) |>
  plot_my_surv2("DLBCL patients OS with MYD88-L265P mutation (n = 1890)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ MYD88_L265P, data = _) |>
  plot_my_surv2("DLBCL patients PFS with MYD88-L265P mutation (n = 730)")

## CD79B gross ----
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ CD79B, data = _) |>
  plot_my_surv2("DLBCL patients OS with CD79B mutation (n = 2153)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B, data = _) |>
  plot_my_surv2("DLBCL patients PFS with CD79B mutation (n = 988)")

## CD79B-Y197H ----
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ CD79B_Y197H, data = _) |>
  plot_my_surv2("DLBCL patients OS with CD79B-Y197H mutation (n = 1890)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B_Y197H, data = _) |>
  plot_my_surv2("DLBCL patients PFS with CD79B-Y197H mutation (n = 730)")

## CD79B Y197S
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ CD79B_Y197S, data = _) |>
  plot_my_surv2("DLBCL patients OS with CD79B-Y197S mutation (n = 1890)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B_Y197S, data = _) |>
  plot_my_surv2("DLBCL patients PFS with CD79B-Y197S mutation (n = 730)")

## CD79B Y197X ---------
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ CD79B_Y197X, data = _) |>
  plot_my_surv2("DLBCL patients OS with CD79B-Y197X mutation (n = 961)")

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B_Y197X, data = _) |>
  plot_my_surv2("DLBCL patients PFS with CD79B-Y197X mutation (n = 730)")

## MYD88 + CD79B ----------
trini_tidy |>
  filter(!CD79B) |>
  survfit(Surv(os_month, os_death) ~ CD79B + MYD88, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             palette = c('blue','orange','green2','red'),
             
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients OS with unspecific MYD88 & CD79B mutation (n = 2153)",
                              width = 55))

trini_tidy |>
  filter(!CD79B) |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B + MYD88, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             palette = c('blue','orange','green2','red'),
             font.legend = 14,
             #legend.labs = c('WT', 'Mut-MYD88', 'Mut-CD79B', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients PFS with unspecific MYD88 & CD79B mutation (n = 988)",
                              width = 55))

## CD79B + MYD88-L265P ----
## 871*590
trini_tidy |>
  filter(!CD79B) |>
  survfit(Surv(os_month, os_death) ~ CD79B + MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             palette = c('blue','orange','green2','red'),
             #legend.labs = c('WT', 'MYD88-L265P', 'Mut-CD79B', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients OS with MYD88-L265P & unspecific CD79B mutation (n = 1890)",
                              width = 55))

trini_tidy |>
  filter(!CD79B) |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B + MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             palette = c('blue','orange','green2','red'),
             #legend.labs = c('WT', 'MYD88-L265P', 'Mut-CD79B', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients PFS with MYD88-L265P & unspecific CD79B mutation (n = 730)",
                              width = 55))

## CD79B-Y197X + MYD88-L265P ------
trini_tidy |>
  survfit(Surv(os_month, os_death) ~ CD79B_Y197X + MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             palette = c('blue','orange','green2','red'),
             legend.labs = c('WT', 'MYD88-L265P', 'CD79B-Y197X', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients OS with MYD88-L265P & CD79B-Y197X mutation (n = 961)",
                              width = 55))

trini_tidy |>
  survfit(Surv(pfs_month, pfs_death) ~ CD79B_Y197X + MYD88_L265P, data = _) |>
  ggsurvplot(pval = TRUE,
             pval.size = 8,
             font.legend = 14,
             palette = c('blue','orange','green2','red'),
             legend.labs = c('WT', 'MYD88-L265P', 'CD79B_Y197X', 'Double-mutation'),
             risk.table = 'nrisk_cumcensor',
             title = str_wrap("DLBCL patients PFS with MYD88-L265P & CD79B-Y197X mutation (n = 730)",
                              width = 55))

# ICGC-GCB -------
icgc.dfs <- read_tsv('mission/DLBCL-Btk/data/icgc-gcb/Disease_free_survival_20240418.tsv')

icgc.dead <- read_tsv('mission/DLBCL-Btk/data/icgc-gcb/donors_2024_04_18_03_38_41.tsv')

icgc.dfs |>
  mutate(group = ifelse(str_detect(donor_set_name, 'no'), 'Not mutated', 'Mutated'),
         dead = donor_id %in% icgc.dead$`Donor ID`) |>
  count(group, dead)

icgc.dfs |>
  mutate(group = ifelse(str_detect(donor_set_name, 'no'), 'Not mutated', 'Mutated'),
         dead = donor_id %in% icgc.dead$`Donor ID`) |>
  survfit(Surv(time, dead) ~ group, data = _) |>
  ggsurvplot(pval = T,
             pval.size = 8,
             font.legend = 14,
             legend.labs = c('Mutated (n=5)', 'Not mutated (n=213)'),
             title = 'GCB lymphona patients DFS with MARCKS mutation')

icgc.os <- read_tsv('mission/DLBCL-Btk/data/icgc-gcb/Overall_survival_20240418.tsv')

icgc.os |>
  mutate(group = ifelse(str_detect(donor_set_name, 'no'), 'Not mutated', 'Mutated'),
         dead = donor_id %in% icgc.dead$`Donor ID`) |>
  count(group, dead)

icgc.os |>
  mutate(group = ifelse(str_detect(donor_set_name, 'no'), 'Not mutated', 'Mutated'),
         dead = donor_id %in% icgc.dead$`Donor ID`) |>
  survfit(Surv(time, dead) ~ group, data = _) |>
  ggsurvplot(pval = T,
             pval.size = 8,
             font.legend = 14,
             legend.labs = c('Mutated (n=6)', 'Not mutated (n=248)'),
             title = str_wrap('GCB lymphona patients OS with MARCKS mutation', 45))

# TCGA-DLBC ------
firehose <- read_csv('mission/DLBCL-Btk/results/TCGA-DLBC-survival.csv')

homodel <-
'dlbc_tcga:TCGA-GS-A9TT-01
dlbc_tcga:TCGA-GS-A9TQ-01
dlbc_tcga:TCGA-FF-A7CR-01
dlbc_tcga:TCGA-GR-A4D9-01
dlbc_tcga:TCGA-FF-8042-01' |>
  read_delim(col_names = c('proj','sample')) |>
  pull(sample) |>
  str_remove('-01$')

firehose |>
  mutate(deleted = patient %in% homodel) |>
  ggplot(aes(deleted, MARCKS)) +
  stat_summary() +
  geom_jitter(width = .1, height = 0)

firehose |>
  mutate(deleted = !(patient %in% homodel)) |>
  survfit(Surv(days, vital_status == 'Dead') ~ deleted, data = _) |>
  ggsurvplot(pval = T,
             pval.size = 8,
             font.legend = 14,
             legend.labs = c('Deleted (n=5)', 'Not mutated (n=43)'),
             xlab = 'Days',
             title = str_wrap('DLBCL patients OS with MARCKS copy number variation', 45))
